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Abstract 

The massive Gross-Neveu model is solved in the large N limit at finite tem- 
perature and chemical potential. The scalar potential is given in terms of Jacobi 
elliptic functions. It contains three parameters which are determined by transcen- 
dental equations. Self-consistency of the scalar potential is proved. The phase 
diagram for non-zero bare quark mass is found to contain a kink-antikink crystal 
phase as well as a massive fermion gas phase featuring a cross-over from light to 
heavy effective fermion mass. For zero bare quark mass we recover the three known 
phases kink-antikink crystal, massless fermion gas, and massive fermion gas. All 
phase transitions are shown to be of second order. Equations for the phase bound- 
aries are given and solved numerically. Implications on condensed matter physics 
are indicated where our results generalize the bipolaron lattice in non-degenerate 
conducting polymers to finite temperature. 



1 Introduction and motivation 



In its original form, the Gross-Neveu (GN) model pQ is a relativistic, renormalizable 
quantum field theory of N species of self- interacting fermions in 1+1 dimensions with 
Lagrangian 

N 1 / N \ 2 

£ = E^(i7^-mo^ (,) + ^ EW 1 • (1) 

8=1 1 \i=l ) 

The common bare mass term oc m explicitly breaks the discrete chiral symmetry ip — > 
^ip of the massless model. As far as the phase diagram is concerned, the 't Hooft 
limit N — > oo, g 2 oc 1/N is most interesting since it allows one to bypass some of the 
limitations of low dimensions, while justifying a semi-classical approach. Interest from 
the particle physics side stems from the fact that the simple Lagrangian (0) shares non- 
trivial properties with quantum chromodynamics (QCD), notably asymptotic freedom, 
dimensional transmutation, meson and baryon bound states, chiral symmetry breaking in 
the vacuum as well as its restoration at high temperature and density (for a pedagogical 
review, see [2]). It also has played some role as a testing ground for fermion algorithms on 
the lattice j^j. Perhaps even more surprising and less widely appreciated is the fact that 
GN type models have enjoyed considerable success in describing a variety of quasi-one- 
dimensional condensed matter systems, ranging from the Peierls-Frohlich model 4j over 
ferromagnetic superconductors [5 to conducting polymers, e.g. doped trans-polyacetylene 
By way of example, the kink and kink-antikink baryons first derived in field theory 
in the chiral limit (mo = 0) [7| have been important for understanding the role of solitons 
and polarons in electrical conductivity properties of doped polymers. Likewise, one can 
show that baryons in the massive GN model (m 7^ 0) are closely related to polarons and 
bipolarons in polymers with non-degenerate ground states, e.g. ds-polyacetylene [HI HUGH] ■ 
Finally, GN models with finite N have been found useful for describing electrons in carbon 
nanotubes and fullerenes [TTJ H2] . 

The claim that a relativistic theory is relevant for condensed matter physics is at first 
sight very provocative [T3]. A closer inspection shows that a continuum approximation 
to a discrete system, a (nearly) half-filled band and a linearized dispersion relation of the 
electrons at the Fermi surface are the crucial ingredients leading to a Dirac-type theory. 
The Fermi velocity plays the role of the velocity of light and the band width the role of 
the ultra-violet (UV) cutoff. In some cases the correspondence is literal so that results 
can be taken over from one field into the other. This was illustrated in ^3] where we 
borrowed results from the theory of non-degenerate conducting polymers, in particular 
the bipolaron lattice, for solving the zero temperature limit of the massive GN model. 

In the process of renormalization the bare mass mo of the GN model is replaced by a 
physical constant 7 > 0, whereas the relation between the coupling constant and the 
cutoff merely sets the scale of the theory. As a result the massive GN model can be 
considered as a one parameter family of physical theories, the point 7 = corresponding 
to the massless (chiral) case. In this paper we determine the full phase diagram of the 
massive GN model as a function of temperature T and chemical potential \x for any value 
of 7. This issue was addressed by Barducci et al. in 1995 ^3] (for earlier, partial results, 
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see also [16J). Recent findings about the massless GN model [13 [TBI US] cast doubts 
on the full validity of their calculations. The assumption that the scalar condensate is 
spatially homogeneous is apparently too restrictive and misses important physics related 
to the existence of kink-antikink baryons. As we now know, the phase diagram in Ref. 
has neither the correct limit m —> 1!) . nor the correct limit T — > |14j . In view of 
the simplicity and generic character of the GN model, its potential for applications to 
different fields of physics, and the fact that it was conceived more than three decades ago, 
it seems desirable to resolve these contradictions. 

The key to the solution of the problem is twofold. First we conjecture the expression 
for the self-consistent scalar potential. To this end a three parameter ansatz will be 
presented and discussed in Sect. 2. Second we formulate the grand potential in a way 
suitable for analytical calculations. In Sect. 3 we show that the grand potential fits into 
the structure of complex analysis. Minimizing it leads to three transcendental relations 
connecting the parameters of our ansatz to the physical quantities T, /i, and 7. From these 
equations it will be possible to read off self-consistency. We outline an efficient procedure 
to numerically solve the parameter equations and present the full phase diagram in Sect. 
4. A detailed comparison with previous results will be given. We prove that all phase 
transitions are second order and investigate singular points on the phase boundary. A 
Ginzburg-Landau effective theory valid in the vicinity of the multicritical point is deduced 
and finally the zero temperature limit will be discussed. Sect. 5 contains a brief summary 
and outlook. In Appendix A we collect some technicalities needed for the proof of self- 
consistency, whereas Appendix B contains the basic formulae necessary for numerical 
calculations. A short version of this paper without detailed formalism has been presented 
previously |20j- 

2 Generalized scalar potential for the massive Gross- 
Neveu model 

A convenient framework for studying the thermodynamics of the GN model is provided 
by the relativistic Hartree-Fock approximation, expected to become rigorous in the large- 
N limit. Due to the (?/^) 2 -interaction in the Lagrangian ([TJ, the Dirac-Hartree-Fock 
equation assumes the form 



with real scalar potential S(x). The scalar potential is related to the thermal expectation 
value of ijjip [see Eq. (J5E]) ]. It thus depends on the eigenfunctions ip and eigenvalues E 
of (J2J) through a self-consistency relation. At present S(x) is regarded as an external 
potential acting on the fermions; self-consistency will be proved in Sect. 3. We choose the 
following representation of the 7-matrices, 






7 1 = ia 3 , 7 5 = 7°7 1 = -a 2 . 



(3) 
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In this representation, the equations for the upper and lower components <f>± of the Dirac 
spinor ip can be decoupled by squaring Eq. ((2)), 

Eq. (j3J) states that the Schrodinger-type Hamiltonians with potentials U± = S 2 ± 5" 
have the same spectra, a textbook example of supersymmetric quantum mechanics with 
superpotential S. Guided by all the special cases already known, we shall use the following 
ansatz for S(x) depending on three real parameters A, k, b 

S(x)=AS(0, Z = Ax, 

5(0 = K 2 sn bsn£sn(Z + b) + Cnb ^ 6 . (5) 

sn b 

Here sn, cn, and dn are Jacobi elliptic functions of modulus k which will be suppressed 
throughout the paper. Since this ansatz is central for the following, we give a detailed 
discussion of its properties as well as its physical origin. 

1) Symmetries 

As a real function, S has a discrete translational symmetry in £ and b with period 2K 
(in this paper K, E denote complete elliptic integrals of first and second kind, K', E' the 
complementary ones with argument k' = \Jl — k 2 ), 

S(t b) = ~S{i + 2K, b) = S(t b + 2K). (6) 

Moreover, it is symmetric under reflections about £ = —6/2 and about £ = K — 6/2 as 
well as antisymmetric under simultaneous reflections of £ and b, 

S(t,b)=S{-b-t,b) = -S{-t,-b). (7) 

From the Dirac-Hartree-Fock equation (J2J) we see that S(x) and —S(—x) are equivalent. 
Therefore we find that the space of different configurations is parameterized by 

A>0, 0<b<K, 0<k<1. (8) 

In the complex £-plane S has a second period of 2iK', 

S(£,6)=£(£ + 2iK',6). (9) 

This, together with the fact that it is meromorphic, makes S an elliptic function in £. (It 
is elliptic in b, too, but we will not need this here.) As an elliptic function it becomes 
worthwhile to study S in the complex plane. 

2) Special values 
From 

S>(0 = /t 2 sn6sn(2£ + 6) 



K 



! sn 2 £sn 2 (£ + 6) 



(10) 
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we see that on the real line S has its extrema at the symmetric points — o/2 + nK; minima 
for even n and maxima for odd n. In the complex £-plane S is of degree two: It has two 
poles (or one double pole) in its fundamental domain. We find 

3{-b/2) = S min = — 

5(K/2-6/2) = 5(3K/2-6/2) = dnb - K ' s f b 

snocno 

5fv h/o\ 5 cno-dno+ 1 

5(K-6/2) = S max = 

sno 

S(iK') = S(-b + iK') = oo. (11) 
The Laurent series at £ = iK' starts with the coefficients 

1 1 + k 2 

c_i = res^iK'S = 1, c = 0, c x = — — . (12) 

sm o 3 



The Laurent coefficients c' k at £ = — 6 + iK' are determined by reflection symmetry 
4 = (-l) fe c fc . 

3) Alternative expressions 



The expression (JSJ for 5 can be cast into various alternative forms. Any of these ex- 
pressions can be easily checked by comparing its pole structure with Eq. (fT2|) . The 
boundedness theorem of complex analysis guarantees that two elliptic functions of equal 
periods are identical up to an additive constant once their pole structure is the same. 

The symmetrized S(^ — b/2) can be expressed in terms of the Weierstrass elliptic function 
V{z) with periods u x = 2iK', u 2 = 2K (see HQ, £± = 6 ± b/2) } 



2 P(£ + + iK') - P(£_ + iK') 

V'{b/2 + iK') P"(6/2 + iK') 



(13) 



V(0 - V(b/2 + iK') 2V'(b/2 + iK' 
or as ratio of Jacobi elliptic functions 



~ , = sn £ + cn £ + dn £ + + sn £_ cn £_ dn £_ 
sn 2 £ + — sn 2 £_ 



Other alternatives for are 



pc / j.\ /t 2 sno cno 

5 (0 = -- r^-cn£cn(£ + 6) + 



Moreover, note that 



dno ' sno dno 

sn , y , ,\ dno 

-dnf dn(f + 6) + — -. 15 

cn o snocno 



^ ( ° 2 = ^b + 1 " ^ " ^ " + b) - (16) 



4) Limits 

The elliptic modulus k varies between and 1. At the boundaries of this interval, we find 

lim,S(£) = cot 6, 

K—>0 

limtf(f) = coth6 + tanh£-tanh(£ + 6). (17) 

The first limit is a constant, the second represents a single baryon profile. The limits 
k, — > and k — > 1 will later describe the transition from the crystal phase with periodic S 
to the massive Fermi gas phase with constant 5* and mass A cot 6 or A coth b, respectively. 
Our ansatz for S is capable of parameterizing both a periodic crystal and a spatially 
homogeneous phase. 

The limiting values for b are 6 = and b = K where we find 

S(£) = 1/6 + 0(6) and 
limS(0 = * 2 ^1 (18) 

Zero is the only value of 6 where k = and k = 1 can coexist and will later be a singular 
point in the phase diagram. The limit 6 — > K yields the form of the self-consistent 
potential familiar from the chiral limit ^H] ■ 

5) Hamiltonian potentials 

From the singularity structure of S, Eq. ()12j) . we read off that the Hamiltonian potentials 
U±(x) = A 2 U±(^), U± = S 2 ± S' have only one pole which is of second order. For U- the 
singularity at £ = iK 7 remains and we obtain the Laurent coefficients 

1 1 + K 2 

c_ 2 = 2, c_i = res ?=iK '^- = 0, c = — g-r — . (19) 

sn z 6 3 

Moreover £/_(£) is an elliptic function and shares periods and poles with 2K 2 sn 2 £. Com- 
parison with Co in the above expansion immediately yields 

t/L(0 = 5(0 2 - £(£)' = 2K 2 sn 2 £ + " 1 " « 2 - (20) 

sm 6 

The result for U + is obtained by reflection symmetry £/+(£) = £/_(— 6 — £). 
#j Spatial averages 

We need the spatial averages of 5 and S 12 in the following sections. While the calculation 
of the latter is straightforward using either Eq. (|16|) or Eqs. (J421 I^Uj) . the first one is 
more intricate. To derive it we have used the Fourier decomposition of the Jacobi elliptic 
functions to get 

1 /- 2K „ 7X 2tt 2 ~ g 2n+1 (2n + l)7r6 , , 



n=0 
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where q = exp(— 7rK'/K) is the elliptic nome of k. Upon comparing the Fourier expansions 
we see that the above sum equals (note a sign misprint in formula 16.23.10 of [23 J) 

f l , d K K2z 

Jfe/K dg sn (zK) 27r^sn o 

with Jacobi's Zeta function Z = Z(b, n) (in the following we suppress the arguments if 
they are b, k). We have finally arrived at 

2K . „, , cnfednfe 



(S) = / d^5(0 = Z 



o sn b 

/•2K 1 TP, 

<S 2 > = / d£ S(0 2 = — h + 1 - - 2-. (23) 

For later convenience we introduce three basic functions of b and k, 

1 cnbdnb E . . 

t = « = 1 - 17' (24) 



sn 2 6' sn 3 6 ' K 

where s and t are related by the elliptic equation 

t 2 = s (s-l)(s- k 2 ). (25) 

The spatial averages can then be written in the compact notation 

(S) = Z + t/s, 

(S 2 ) = s-1-k 2 + 2u. (26) 

Averages over higher powers or derivatives of S can also be evaluated in closed form, see 
Eq. flUD. 

7) Physical origin 

We close this section with a short digression on the role this form of the scalar potential 
plays in condensed matter physics and how it was deduced in this context. We first became 
aware of this ansatz through works on the bipolaron lattice for non-degenerate conducting 
polymers j2I], [22] , |2*l]-|27j. The following discussion is based on these references as well 
as on mathematical literature on periodic solitons |2~H1 12"§] . 

Let us first recall the results for single baryons in the massive GN model [HI Elj (or, 
equivalently, polarons and bipolarons in quasi-one-dimensional condensed matter systems 
jHj). The scalar potential for a baryon has the form 

S(x) — 1 + y [tanh(yx — Co) — tanh(yx + Co)] , Co = - artanhy. (27) 

The parameter y depends on the bare fermion mass and the number of valence fermions. 
Here the corresponding isospectral potentials U± in the second order equation (JJJ) are 
given by the simplest Poschl- Teller potential [30J 

g2 _|_ g> ^ ^28) 

cosh 2 (yx ± Co) 
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and differ oniy by a translation in space (they are "self-isospectral" in the terminology of 
Ref. [SI]). Their distinguished feature is the fact that they are reflectionless; indeed, they 
are the unique reflectionless potential with a single bound state. It is well known that 
static solutions of the GN model must correspond to reflectionless Schrodinger potentials 
E2] . The fact that the ansatz (}2*T|) leads to self-consistency as shown in Refs. [TUJ is 
therefore quite plausible. 

Take now a lattice of infinitely many, equidistant Poschl- Teller potential wells. As dis- 
cussed in Refs. [221221, the lattice sum can be performed yielding a Lame-type potential, 



E 



i 



'2kK!' 



2 r 



cosh [x — nd) 



E' 



/2K' 



sn 



-x 



7T 



Comparing the spatial period of both sides of Eq. (1291). we can relate d and k via 



d = 7T: 



K 
K 7 



(29) 



(30) 



How does the fact that the single potential wells (|28J) are reflectionless manifest itself in 
the periodic extension (|29*j) ? This has been discussed in mathematical physics [2E] and 
condensed matter physics (32] some time ago: The periodic potential has a single gap 
(or, in general, a finite number of gaps), in contrast to generic periodic potentials with 
infinitely many gaps. Thus reflectionless potentials generalize to "finite band potentials" 
as one proceeds from a single well to a periodic array. In the same way as the sech 2 - 
potential is the unique reflectionless potential with one bound state, the sn 2 -potential is 
the unique single band potential. Guided by these considerations, let us try to find the 
most general superpotential of the Lame potential (plus an additive constant). Allowing 
for a translation in space between U± and using an appropriately rescaled coordinate £, 
we have to solve the equations 



U+ = S 2 + S' 
u_ = S 2 -S' 



2/€ 2 sn 2 (£ + 6) +r], 
2/€ 2 sn 2 £ + r] 



(31) 



for S or, equivalently, [compare Eqs. (fTUl [TBjl ] 

S 2 = K 2 [sn 2 (£ + 6) + sn 2 £]+r7, 
§' = K 2 [sn 2 (e + 6) -sn 2 £]. 

Differentiating the upper equation using 

(snf)' = cn£dn£ 



(32) 



(33) 



and dividing the result by the lower equation yields exactly the ansatz (jSJ) in the equivalent 
form (JHJ). As a by-product, by specializing Eq. (|52Jl to £ = we can determine the 
constant r] [cf. Eq. (|2T!J) ]. 

1 



sn 2 6 



1 



(34) 



The scale factor A in Eq. (jSj) is not constrained by these considerations. Thus we conclude 
that the ansatz (jSJ) is the most general Dirac potential leading to a single gap Lame 
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potential (plus constant) in the corresponding second order equations. This makes it a 
good starting point for finding periodic, static solutions. 

It could be shown in ^1] that the ansatz (j!4)) describes the zero temperature limit of the 
massive GN model. Regarding S(x) the only novelty is now that all three parameters A, 
K, b will be varied independently whereas at T = the parameter A was given by the 
density and only k, b were to be determined by a minimizing procedure. In the chiral 
limit m = we already found in that the form of the scalar potential at T = 
extends to finite temperature without any other change. This result encouraged us to 
follow a similar strategy in the massive case. 

The only way of further generalizing the ansatz would be to go to higher genus Lame 
potentials, increasing the number of bands, but for physical reasons we do not expect 
this to be necessary to describe the thermodynamics of the GN model. Notice that 
the simple and intuitive relation between the single baryon and the crystal exhibited in 
Eq. is hidden in the corresponding Dirac potentials (|27)l and © due to the non-linear 
relationship between S and U±. 



3 Minimization of the grand canonical potential and 
proof of self-consistency 

In this section, we employ the ansatz ((5J) for the scalar potential S(x), minimize the grand 
potential with respect to the three parameters A, b, k and prove the self-consistency of 
the scalar potential. 

Collecting the information derived in the previous section we have transformed the Schro- 
dinger type eigenvalue equation (jlj) into the single gap Lame equation 

/ d 2 \ 

+ 2/t 2 sn 2 £ 0+ = £<f) + , (35) 



and a shifted (£ — > £ + b) equation for 0_. The Lame eigenvalues £ are related to the 
Dirac energies E through 

e = 5 - \ + 1 + k2 - ( 3e ) 

A z smo 

The solutions of Eq. (J35|) are well known, see e.g. Refs. jHHEH]- For logical completeness we 
will sketch the results to the extent needed and follow the method used in the chiral limit 
PHI to transform the grand potential and the thermal expectation value of ipip into complex 
contour integrals. Already in the previous section it turned out to be advantageous to 
extend functions into the complex plane and use results from complex function theory. 
Whereas in JH] such an extension was used as a welcome simplification, here this step 
is vital to solve a much more intricate problem. It is remarkable how well the physically 
real problem fits into the language of complex analysis. 
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1) Eigenfunctions and eigenvalues 

We review the properties of the normalized (dimensionless) eigenfunctions of Eq. (|35|) 
which by abuse of notation are again denoted by ip = (<f>+, <P-) T ■ 

It is possible to give a representation of + in terms of H, 0, and Z, the Jacobi eta, theta, 
and zeta function [T9"] . 

M0=^^^-^ t \Af\ 2 = K f (37) 

B(£) 7rB^(a)|dn a — E/K| 

To any Lame-eigenvalue £ there exist two complex conjugate eigenfunctions parameterized 
by ±a. For ±a = K + iK' ... K one obtains the solutions of the lower band k 2 < £ < 1 
whereas ±a = ... iK' provides the upper band £ > 1 + k 2 . 

The eigenfunctions are presented in the form anti-periodic function times complex phase. 
This gives them a definite Bloch momentum which is opposite for a and —a. The normal- 
ization M is chosen such that the spatial average obeys (|0+| 2 ) = 1/2 [see Eq. (J32J)]. The 
energy eigenvalue together with the Bloch momentum give a parametric representation 
of the dispersion relation, 

£ = 1 + K 2 cn 2 a, p = -iZ(ex) ± (38) 

where we had to add the offset ±7r/(2K) to the momentum to give the periodic lowest 
energy solution a = ±(K + iK') zero Bloch momentum. From the dispersion relation we 
can calculate the density of states and determine its high momentum asymptotics, 

dp \£-k 2 -E/K\ 




d£ 2y]{8 -k 2 )(£ -1){£-1-k 2 ) 

£ = p 2 + 2(l-E/K) + 0(p- 2 ). (39) 

It is worthwhile noting that S(£) = — ln0 + | a= _ fe 
With cf) + given, 0_ follows from the Dirac equation (J2J), 

(40) 

From this equation we see that 0_ has the same Bloch momentum as 0+ (as physically 
demanded). Moreover we know that 0_ solves the 6-shifted Lame equation. This makes 
</>-(£) a multiple of 0+(£ + b). To determine the ratio we observe that the reflected and 
charge conjugated Dirac spinor ■?/;(£) = — 7oV ; ( — b — £)* = (0^— 6 — £), 0^—6 — £)) T fulfills 
the Dirac equation (j2J for the same energy eigenvalue as ip and is thus a linear combination 
of ip and ifj*. Because 4>±(— b — £)* have the same Bloch momentum as 0±(£) we conclude 
that ip is in fact proportional to ip. In components this means <p*_(—b — £) = c0 + (£) and 
0*j_(— 5 — £) = c</>_(£). Upon iteration we find cc* = 1. From the explicit form of + , 
Eq. (|3*7J). we know that 0+(— £)* = ±</>+(£), depending on the energy band. We conclude 
that </>-(£) differs from 0+(£ + b) by a complex phase, 

0_(£) =6^0+ (£ + &), with^GM. (41) 
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This makes the spatial average of their squared absolute values equal and we obtain 

^V) = (|0 + | 2 ) + (|0-| 2 )=2(|0 + | 2 ) = 1 (42) 

as demanded. 

Finally we calculate ^0(0 an d V^VKOj the sca l ar an d baryon densities of a single orbit. 
We find |E| 

^=>* +25 >^ 2 > i^r= L £ :;::g - <«> 

We express £ in terms of E (Eq. EEJ), use the identity 

I , 9j _ sn 6 cn £ dn £ + sn £ cn 6 dn b ,,,, 
- 1 + dn 2 £ = 9 , s 77 ^V. (44) 



sn 2 6 H sn 2 6sn(6 + £) 

and Eqs. (J2U) to arrive at the result 



= (^)j-v(^) . (45) 

[El Ay — s + u 

From (gU) we know that f 0(£) = 10+ (01 2 + l<Mf + h )? ■ With Eqs. (HEJ) and (USD this 
yields 

^ Y (E/A) 2 -s + u K ' 

2) Grand canonical potential, scalar condensate, and baryon density 

In Hartree-Fock approximation, the grand canonical potential density per flavor consists 
of the independent particle contribution \l>i and the double counting correction \1> 2 , 

*a = --r- / A/ dgln[(l+e-^-^)fl + ^+^ 
pn Jo LV ' v ' 

^2 = TTTT^ f £ MS(x) - m f. (47) 



2NgH Jo 

Here, A/2 is the UV cutoff, £ = 2K./A denotes the spatial period of S(x) and q is the 
Bloch momentum. The fact that mo 7^ manifests itself in two places. Via S(x) it 
influences the single particle energies E in and it leads to the shift S(x) — > S(x) — itlq 
in ^2- Let us first transform the momentum integral in into an integral over Dirac 
energies. Scaling out the factor A via 

P-l, «-f (48) 
and using Eqs. (|24j) the change of integration measure and cutoff is given by [see Eqs. (J3bl 

mi 

dp UJ(UJ 2 -S + U) . 2 1U 2 2 \/ 2 \ 

dcu ±\/W 

A. = w (A/2) = ^ + ^ + °( A_3 )- ( 49 ) 

ll 



The plus sign in front of v W refers to the upper band, the minus sign to the lower band. 
Due to the quadratic divergence of the integral, it is necessary to expand A w to order A -1 . 
With the shorthand notation 

a = (3A, u = (3n (50) 



can be written as the following integral over the energy bands \J s — 1 < uj < Vs — k? 
and uj > y/s, 



vr/? 2 ^! = -a ( [ * *duj + f^duj) ^ In [(1 + e~ auj+u )(l + e auJ+u )] . (51) 

Following PHI , we combine the integral over both energy bands (including the sign change 
in VW) as well as over positive and negative energy modes to the real part of a line integral 
in the complex cu-plane. The path of integration runs above (or below) the real axis and 
the sheet of yW is chosen such that \/W\ u=+o0 = +oo (see Fig. 1) yielding 

tt/? 2 ^! = -aRe lim du ^ ~* + U) In f 1 + e^ +") . (52) 



+ 1 



K 



-\/s - I \fs 



Figure 1: Path of integration in the complex u;-plane running slightly above the real axis. 
The fat lines indicate branch cuts of \/W in Eqs. (I49| I52"J) corresponding to the gaps in 
the Dirac-Hartree-Fock equation We have indicated the complex phase of \fW above 
the real line. Note that \/W is antisymmetric under reflections about the imaginary axis. 



The contribution of ^>2 was calculated in 

* 2 = ((S 2 ) - 2m (S) + ml) 

= ^^{ a2 ( s - 1 - K2 + 2u )- 2a ^ m o^ + t/s)+(3 2 m 2 ). (53) 

With the technique used to rewrite \l/ we transform the thermal expectation values of the 
scalar potential and the baryon density into complex integrals. Using Eqs. (|45|) and (|49J) 
we obtain 

1 r A / 2 - ( 1 1 \ 

W) * = nJo dg ^ U^) + l ~ e-^) + J 

a _ .. r°°+ ie , uj 2 S - 1 1 , rA , 

= — Rehm / do; = , (54) 

7T/3 e^oJ-W +i£ y/w e^-^ + l v ' 
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and analogously with Eq. (|4*H|) . 

p(x) = ma = A R e lim ^ + (g 3 -3 f _+l + « 3 )/2] 1 , (55) 



Because the integrals in (|54|) and f)55|) are only logarithmically divergent it is sufficient to 
keep the leading term of K w . With Eqs. (|52p and (J53j) we confirm that the spatial average 
of the baryon density satisfies (p(x)) = —d^f/d/j,. 

3) Self- consistency 

We have to minimize \& with respect to A, b, k for fixed /x, /3, 7 and prove that at the 
minimum the self-consistency condition 

holds. 

It is convenient to rescale the integration variable in (J52j) by u — > uj/a before differentia- 
tion. The contour of integration stays off the singularities of the integrand by a minimum 
distance of e (the limit e — > is only needed at the far left end of the integration contour). 
This allows us to interchange differentiation and integration and we obtain equations with 
integrands of the generic form (polynomial in uS) x iy~ 3//2 ln[l + exp(— auj + v)\. Noting 
that <9 w ln[l + exp(— au + u)) = —a/[exp(auj — v) + 1] we see that we gain an equally 
looking expression upon applying a partial integration on the scalar condensate given by 
Eq. ()54|) . Comparing coefficients of the polynomials such that the self-consistency condi- 
tion (}56|) holds on the level of integrands yields an over-determined system of equations 
which happens to have a solution. This solution actually gives the full self-consistency 
condition. 

In order to present this solution in a most concise way we use the freedom to minimize 
\I/ with respect to any set of independent variables. Specifically, we propose to replace 
A, b, k by the spatial averages (S), (S 2 ) and the spatial period £ of S. The stationarity 
conditions for the grand potential then read 

d 

Tip 2 ^ = (3F = 0, 



d(S) 



9 o2,r, P 



2 



d{S 2 ) r 2 

^ - °jF, - 0, (57) 

where we have introduced functions Fi which play a key role in the following. The factors 
in front of the Fi have been included for later convenience. In the limit tuq — > we find 
that Fx and F2 equal the corresponding functions in ^H] • The derivatives with respect to 
the new, composite variables can be taken trivially in the case of ^2, 

d(S) Ng 2 ' d{S 2 } 2Ng 2 ' 3£ 1 j 
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Unfortunately this is not true for available only in terms of the original variables A, 
b, k in Eq. Q52J). In order to compute Fi, we therefore invoke the chain rule, 




/ d(s)_ d(s 2 ) m \ 

Ah Ah Qh 

M 
c)A 



db db db 

d(S) d(S 2 ) 

dA dA 

\ d(s± d{s^ ai I V a 2 F 2 /i 

\ 8k dn 8k I 



(59) 



Upon inverting the Jacobian matrix on the right hand side, the F{ can be expressed in 
terms of more accessible derivatives of Details are given in Appendix A. We arrive at 



F 
F-2 



7r/3m 

"IV" 



/oo+ie I Q 

da; 



7T 



—Re lim 

a e^O 



oo+ie 
oo+ie 



Ng< 



I roc+ie 

—Re lim / duj 



did yfw ) 
< 8 u 2 



In 1 



0. 



dw _ 

+ie V # 



In 1 



0. 



tufa; — s + w) 



d tZ- ULO 2 + (lo 2 - s + 1)(cj 2 - S + K 2 



duo 



ln(l + 



-auj+u 



0. 



(60) 



Note that the integrand in F 2 falls off at u — > ±oo fast enough to render the integral 
convergent. All three functions Fj vanish at the minimum of the thermodynamic potential 
(|52 | I53 |) . After a partial integration we compare F and F with Eq. (J54)) and obtain 



(61) 



At the minimum of the grand potential we obtain self-consistency f!56jl from ^0 = ^1 = 0. 
We do not need the minimum condition F 2 = to ensure self-consistency, nor do we 
need the limit In A — >• 00 as long as we drop terms that are 1/A suppressed. This latter 
property becomes significant in condensed matter physics where A is the finite band width 
of the system under consideration. 



4 Phase boundaries and phase diagram 

The main result of this section is the phase diagram of the massive GN model. The dia- 
gram is plotted in Fig. 3 and analyzed in Subsect. 4.4. Before that we have to renormalize 
the model and solve the minimum equations Fi = 0. We end this section with an analysis 
of the phase boundaries and its singular points. 

1) Renormalization 

In the version presented so far the massive GN model still contains an UV cutoff A. For 
a quantum field theory this state is not acceptable, it is necessary to remove the cutoff 
and present a renormalized version of the model. In the wake of this procedure the bare 
parameters Ng 2 , mo, and the cutoff A are replaced by physical quantities that have to 
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be adjusted to observables. Since all UV divergences are due to vacuum effects, there are 
no new difficulties as compared to the T = case [Hj. All we need is the vacuum gap 
equation 

J_ = i(l +mo )l nA= I(l nA + 7 ). (62) 

JM g Z 7T 7T 

Renormalization of Ng 2 interacts with the cutoff A to merely set the scale of the theory. 
We fix this scale by the condition that the dynamical fermion mass in the vacuum equals 
1 . The bare fermion mass mo however is renormalized to a physical parameter 7 = mo In A 
(called "confinement parameter" in condensed matter physics). For different values of 7 
one obtains different physical theories, only one of which would be accessible in a fictitious 
1+1-dimensional world. With 7 > we trace through a one parameter family of physical 
theories where 7 = corresponds to the chiral, massless GN model. Renormalization can 
be performed by sending the cutoff A to infinity while carefully reducing Ng 2 and m such 
that the dynamical fermion mass and mo In A are kept at the values 1 and 7, respectively. 
After renormalization N, g, mo, A have disappeared from the theory in favor of the single 
parameter 7. 

In the grand potential irrelevant divergent terms — A 2 /87r and —fiA/2n stemming from the 
energy and baryon density of the Dirac sea can simply be dropped. A short calculation 
using the antisymmetry of \/W and the fact that due to its holomorphy in the upper 
half plane uo{uo 2 — s + u)/ \/W — 1 vanishes if integrated above the real line shows that an 
expression for the renormalized grand potential is (A = /3A) 



2,T,ren 



lim 

A— >oo 



A 2 a 2 (S 2 ) / A 

1 — - hi - + 7 

8 2 I p ' 



1 



aP 7 (S) 



a Re lim / 2a duo^-^- — ?JL^ U ^ In (2 cosh 

e^o 



(63) 



Note that we have not only symmetrized the integration domain but also replaced A^ by 
its leading term and corrected the effect by subtracting a 2 {S 2 )/2. In fact, the right hand 
side of Eq. (JBlIj) converges with 0(A~ 2 ) as A — > 00. 

The effect of renormalization on Fq and F% amounts to using the gap equation (|62|). 
whereas F 2 is unaffected by renormalization. Due to symmetry (needed for F{) and 
holomorphy in the upper half plane we may replace ln[l + exp(— auo + v)] by the more 
symmetric expression In 2 cosh[(au;— v)/2\. After partially integrating F and F\ we obtain 



F™ = -/3 7 + at/ 



F? 



0. 



7 -ln/3- h = 0, 
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OJ 2 — s 
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In [ 2 cosh 



2 
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(64) 



Note that I , Ii, and F 2 have no explicit /^-dependence. 
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Although we will not make use of it in the following, let us mention that one can use 
Fi = to present \I/ at its minimum in a somewhat simplified way. One obtains 



P mm 8 s 2 



where the limits e — > and A — ► oo are understood. 




W au — v 
tanh 



(65) 



After solving Eqs. (JMj) we can use Eqs. or (J63j) to calculate the pressure P = — ^ as 
a function of p, T giving the equation of state. Bulk thermodynamical variables like the 
(spatially averaged) baryon density p, the entropy density a and the energy density e can 
then be obtained via standard thermodynamic relations, 

P=-q-*> a = ^df 1 ' e = Ta-P + pp. (66) 
2) Solving the minimum equations 

With the physical parameters 7, /3 = 1/T, v = p/T given we have to solve the three 
equations Fq 611 = F± en = F 2 = 0, Eqs. (jMJ), simultaneously to obtain the values of the 
parameters a, b, k at the minimum of the grand potential. There exists a simple and 
efficient method to achieve this numerically. 

We observe that F2 depends on only four of the six variables involved, namely on a, b, 
k, v. For convenience we fix a value of k with < k < 1, of b with < b < K, and 
the physical v. With k and b given, we calculate s, t, u, Z and consider the function 
F 2 (a). We find that F 2 (a) = ca 2 + C(a 4 ) for a -> [see Eq. flUJ)] and that F 2 (a) goes to a 
negative constant for a — ► 00. If (here ip(n, z) = d™ +1 lnT(2;) is the polygamma function) 

v>v t , with Reip (2, - + — ) = 0, v t = 1.910668 (67) 

V 2 2tx J 

we have c > and the function F 2 (a) has a zero at a unique positive value of a. For 
f < u t we have c < and F 2 (a) = has only the trivial solution a = 0. In the latter case 
the crystal solution will not exist and we are in the (simpler) translationally invariant 
regime studied in the following subsection. At the moment we restrict ourselves to the 
case v > Vf 

With F 2 (a) given as a one-dimensional integral over elementary functions the positive zero 
will quickly be found. Now we can calculate the integrals Io and I\ in Eqs. (I64J) which 
also depend on a, b, k, v only. Combining Fq = and F\ = 0, we obtain 

7 + ln7 = ]n(atio)+Ji, P = atl /y. (68) 

These equations determine 7 and 3 for given values of b, k, v. Normally, one would 
like to dial 7, p, T and evaluate thermodynamic observables for these values. Therefore 
one has to work one's way to the desired values of 7 and Q by systematically estimating 
parameters b and n to start with. If 7, p, T lie in the crystal region depicted in Fig. 3 the 
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iteration will terminate quickly. Only at the phase boundary this procedure becomes less 
efficient, but there it will be possible to derive alternative methods by virtue of analytical 
techniques (see Subsect. 4.5). 

We close this subsection with the remark that for actually implementing this algorithm 
on a computer one should transform the contour integrals into real integrals of smooth 
functions with compact support. This is always possible as worked out in detail for a 
variety of integrals in Appendix B. 

3) Homogeneous phase and the "old" phase diagram 

The translationally invariant, homogeneous phase has been studied in |15| . We will briefly 
review the results and fit them into the context of this more comprehensive approach. 

We saw in Eq. (117)1 that in the limits k — > and k — > 1 the scalar potential S tends to 
a constant. The single baryon in the case k — > 1 is insignificant in entire space. We see 
from the Dirac-Hartree-Fock equation (J2J) that in this case the scalar potential has to be 
interpreted as an effective fermion mass M. We have 



S(x) — M — ^ILJl for K = o and M = ^ for « = 1. (69) 

The homogeneous phase is accessible in our framework by setting k = or k = 1 and 
not demanding F 2 = 0. Because the minimum condition F 2 = was not needed for self- 
consistency we will not lose self-consistency in the homogeneous phase. The two bands 
collapse to a single band and the grand potential in either limit specializes to [after a 
rescaling u — > uj/3/a and a partial integration in Eq. 



^hom = lim 

A^oo 



A 2 M 2 1 /■# + ie > 8u-v 

1 (In A + 7) -7M + -Relim / dcuVcu 2 -M 2 tanh 

8 2 2 e^oj-A +; e 2 

(70) 

The limit A — ► 00 can easily be performed. For comparison with other work we give the 
result in the original variables, momentum q and energy E = y/q 2 + M' 2 , 



,J 2 



horn 



(21nM-l) + ^(M-2)- J- I" dq\n\(l + e-^ E ~^) fl + e 
4tt 2ty pir Jo -' ; 



M 2 , 7 Af, , 1 

- / agin lite r> ~ r, i lite ri " ,r ' 

(71) 

The minimum condition is obtained by differentiation of Tr^om with respect to M. Com- 
parison with Eqs. (}64|) [or Eqs. (|T4| I7H|) ] yields 



rcn 



= M In M + 7 (M - 1) + M | ^ ^^^^^ + ^^j^j = -f- + MFJ™. (72) 

In fact, only M, a combination of a and 6, is of relevance. This reduces self-consistency 
to a single linear combination of Fq and Fx. A second equation is not needed. For the 
moment we study the homogeneous phase in the whole (/x, T)-plane, keeping in mind that 
the phase diagram will change due to the existence of the crystal phase. 

Depending on the parameters (7,/i, T), there may be one or two local minima with the 
possibility of a first order phase transition. Since we could not find a published graph 
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1 

1.2 

Figure 2: "Old" phase diagram of the massive Gross-Neveu model as a function of 
(7, //, T), obtained under the (unjustified) assumption of unbroken translational sym- 
metry. The various phase boundaries are explained in the main text. 



of the phase diagram in (7, /1, T)-space, we have recalculated it with the result shown in 
Fig. 2. 

In the chiral limit (7 = 0), the massive phase at low (/x, T) is separated from the chirally 
restored, massless phase at high (/i,T) by a critical line AP t B |36| . The upper part AP t 
of this line is second order, the lower part P t B first order, with a tricritical point P t 
separating the two. As the parameter 7 is switched on, the second order line disappears 
in favor of a cross-over where the fermion mass changes rapidly, but smoothly. The first 
order line on the other hand survives, ending at a critical point. The effect of the cross- 
over shrinks to zero when the critical point is approached. If plotted against 7, these 
critical points lie on the third curve P t C emanating from the tricritical point. For 7 > 
the effective fermion mass M is positive everywhere. If one crosses the shaded critical 
"sheet" in Fig. 2, the mass changes discontinuously, dropping with increasing chemical 
potential. More details on the thermodynamic implications of such a phase diagram can 
be found in Ref. [To] . 

4 ) Crystal phase and the revised phase diagram 

According to the discussion of S(x) in Subsect. 2.4, at the phase boundary between ho- 
mogeneous and crystal phases, k has the value or 1. Only for these limiting values of the 
elliptic modulus, the spatial modulation of S(x) becomes insignificant. This observation 
is the key for deriving the phase boundary equations in the following subsection. We have 
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computed points of constant v — fi/T varying b from to 7r/2 or to oo, respectively. For 
K — or 1, varying v > u t [see Eq. (JH7|) ] we obtain (6, z/)-grids which are mapped onto 
curved two-dimensional surfaces in (7, /1, T)-space. These surfaces represent second order 
phase boundaries (see Subsect. 4.5). 




pppp- 



C 



0.6 

7 0.8 



Figure 3: Revised phase diagram of the massive Gross-Neveu model. The shaded surfaces 
/, 77 separate the kink-antikink crystal from a massive Fermi gas and correspond to second 
order phase transitions. 



We now turn to the results, focusing on the phase boundaries as depicted in Fig. 3. In 
the 7 = plane, the result for the revised phase diagram including the kink-antikink 
crystal is known [T^j: The second order line AP t of the old phase diagram is unaffected. 
The first order line P t B of Fig. 2 is replaced by two second order lines PtB' and PtD 
delimiting the kink-antikink crystal phase. The tricritical point is turned into another 
kind of multicritical point (called Lifshitz- or Leung point in condensed matter physics 
|3T| |38j), located at the same (/i, T) values. As we turn on 7, the second order line 
separating massive (M > 0) and massless (M = 0) phases disappears as a consequence of 
the explicit breaking of chiral symmetry. The crystal phase survives at all values of 7, but 
is confined to decreasing temperatures with increasing 7. For fixed 7, it is bounded by two 
second order lines joining in a cusp. The cusp coincides with the critical point of the old 
phase diagram but has a significantly different character. The crystal phase exists and is 
thermodynamically stable inside the tent-like structure formed out of two sheets denoted 
as / and //. These sheets are defined by k = (/) and k = 1 (//), respectively. The 
line P t C where they join corresponds to b = and coincides with line P t C in Fig. 2. The 
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baseline of sheet II in the (/x, 7)-plane has a simple physical interpretation: It reflects the 
7-dependence of the baryon mass in the massive GN model, or equivalently the critical 
chemical potential at T = 0. The chiral limit 7 — >■ can be identified with b — > K. In this 
way, one can recover the simpler 2-parameter ansatz for the self-consistent potential used 
in from the 3-parameter ansatz (jSJ). Incidentally, we should like to point out that all 
the fat lines in Fig. 3 representing the boundaries of the sheets /, II have been taken from 
sources independent of the present work [THJ HZl HH| • The fact that the computed surfaces 
connect very accurately constitutes a test for our algebraic and numerical calculations. 



0.3- 
0.25- 
0.2 

T 

0.15- 

0.1- 

0.05- 

0.5 0.6 0.7 0.8 0^9 i l!T 1.2 

/' 

Figure 4: Two-dimensional sections 7 = const, through the phase diagram of Fig. 3. The 
fat lines correspond to 7 = 0, the thin lines to 7 = 0.1. ..1.2 (in steps of 0.1), from top to 
bottom. The position of the cusp agrees with the critical point in the old phase diagram. 



Another view of the phase diagram is shown in Fig. 4. Here we plot the phase bound- 
aries in the (/x, T)-plane for several values of 7, i.e., for several bare quark masses. These 
two-dimensional graphs correspond to cutting the three-dimensional graph in Fig. 3 by 
equidistant planes of constant 7. They provide a better view of how the two critical lines 
are joined in the cusp and exhibit that the region in which the crystal is thermodynami- 
cally stable shrinks rapidly with increasing 7. 

In order to further illustrate the nature of the phase transitions between inhomogeneous 
and homogeneous phases, let us remark that outside the tent the scalar potential is x- 
independent, i.e., an effective mass. If one crosses the tent from side / to side // on some 
curve, the parameter k varies continuously from to 1. By way of example consider the 
straight line 7 = 0.5, T = 0.05/i which pierces sheet I at fi — 1.1 and sheet // at /i = 0.85. 




20 



Figure 5: Evolution of self-consistent S(x) as one crosses the crystal region on a straight 
path. As /j, decreases from 1.1 to 0.85, the elliptic modulus k varies between and 1. 
According to the "old" phase diagram, there would be a discontinuous jump between two 
constant values of S at fj, — 0.87. 



Just outside of the crystal phase, the effective mass values we find are 0.40 (7) and 1.0 
(77). According to the old phase diagram, the mass jumps suddenly from 0.54 to 1.0 
upon crossing the first order sheet at /i — 0.87. By contrast, Fig. 5 shows the continuous 
evolution of S(x) in the exact calculation. At k — 0, an instability occurs with respect 
to oscillations of a finite wave number. As k increases, the amplitude of these oscillations 
become larger whereas the period first grows rather modestly. When k approaches the 
value 1 the period grows rapidly and, save for distant dents representing widely spaced 
baryons, S(x) reaches a constant value connecting to the translationally invariant phase. 
At k — 1 the system is instable against single baryon formation. This subtle interpolation 
of S(x) between two constants caused by the Peierls instability is only crudely modelled 
by a translationally invariant scenario. 

Before we examine the phase boundaries and its special points in the next three subsec- 
tions, let us summarize the thermodynamics of the massive GN model in Table 1. It is 
important to notice that the sheet of first order phase transitions present in the old phase 
diagram is entirely superseded by the crystal phase. The grand potential of the GN model 
is smooth for any values of 7 (including the value zero) everywhere in the (/i, T)-plane. 
Even approaching the critical points from any direction does not result in a kink in the 
grand potential: The effect of the cross-over as well as of crystal formation reduces to 
zero when the critical point is reached. 
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Eqs. (mm 
Eq S . Annum) 



Surfaces 



Eqs. fUEl, 
Eqs. (plITgl). 



2 nd order phase boundary 
2 nd order phase boundary 



Lines 



7 = 0, ln(47rT) + ReV(| + f) = 

7 = 0, f I£ /2 ^Imln T M + ^ +aCOS<P) l = 0, 

I ' daaJU cosip r( f + ^r iy-a cos tp)j 

ln(47rT) + i # Re ^ (| + £ + a cos ^)) = 
7 = 0, ln(47rT) + \ minRe Q + + a)) + ^ (§ + i< 
T = 0, + 7tan# = tt/2, sin # + 7 artanh(sin 9) = fiir/2 
Eq. (JEHD 



z/ — a 



Points 



(7, //, T) = (0, 0, e c /vr) = (0, 0, 0.566932), C = 0.577215 (Euler constant) 
(j,fi,T) = (0,z/ t T t ,T t ) = (0,0.608221,0.318329), 
Eq. (EZD, ln(4vrT t ) + Re^ (§ + = 
(7,/x,T) = (0,2/tt,0) = (0,0.636619,0) 



Table 1: Guide through the phase diagram of the massive GN model, using results from 
PI EH HI]. m this table T(z), i[>(z), and ip(n, z) are the Gamma, digamma, and polygamma 
functions and v = [x/T. 



5) Phase boundaries 

We saw in Eq. f!17|) that we lose the crystal shape of S(x) in the limits k — > and 
k — ► 1. In the first case the amplitude of the modulations vanishes, in the second case 
the frequency of the modulations goes to zero. Corroborated by numerical results and 
the 7 = limit in [19] we expect to reach the phase boundaries when k tends towards 
its extremal values and 1 (the boundaries are labelled I and II, respectively, in Fig. 3). 
More precisely, only in the limit 7 = we will have two phase boundaries requiring at 
least three different phases. For 7 > we will find that the only phase boundary between 
the crystal and the massive Fermi gas consists of two different sections separated by a 
distinguished point. This point, the cusp, will be examined in the next subsection. Here 
we focus on the behavior of the theory for values of k close to or 1. 

In the limit k — > the two energy bands combine to a single band uo 2 > s — 1. A remnant 
of the shrunken gap is a singularity at uj 2 = s inside the energy band. It is useful to 
change variables from b to s = sn~ 2 b which gives the expansions in k a compact form, 

8-1 , (5-1X5-2) | (3-1)^-2) ^ 

tL ~ ~^2~ K + lfo + 32^ 
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(73) 



Expanding Fj yields [see Eqs. (JMJ)] 
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(74) 



Solving these equations following the procedure of Subsect. 4.2 with k = leads to a 
relation between the physical parameters 7, /1, T which determines phase boundary /. 

In the limit k — > 1 the procedure is analogous. The lower energy band collapses to a 
point singularity at u 2 = s — 1. We are left with a single band co> 2 > s with no singularity 
inside. Expanding in k' = yl — k, 2 yields 
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Expanding Fj at k = 1 yields 
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These equations, too, can be solved with the procedure of Subsect. 4.2 yielding the relation 
that defines phase boundary //. Neither in the case k = nor in the case k — 1 it is 
possible to formulate the phase boundary condition in terms of the effective fermion mass 
M [Eq. (J69|) ] alone. For numerically accessible versions of the integrals in Eqs. (J741 17U|) 
see Appendix B, Eqs. (ITMITUoT) . 

In the vicinity of the phase boundaries I and II, one can evaluate the difference of \l/ ren 
in the crystal and in the homogeneous phase analytically. The results are 

IT-Tol 



^rcn 

cryst 



j ren 
horn 



-C7(T-T 



cryst 



j ren 
horn 



c 



11 



In IT - T n 



(77) 



23 



respectively, where To is the temperature at which the phase boundary is crossed. In 
both cases, Cj and C/j are positive constants leading to *&™ st < ^hom- The crystal is 
thermodynamically stable near the phase boundary, the phase transition is second order. 

In order to prove Eq. (J77|) and to determine the constants one has to specify the direction in 
which one enters the crystal phase. For simplicity we assume that we stay on a v = const, 
line when passing through the phase boundary. In this case we can use expansions (|75jl 
and (jZSj) to calculate the K-dependence of a, s or b, and f3 from Eqs. (|51jl. Although at 
k = one best works with s instead of b we use s = 1 — sn~ 2 (ife, k') and the last identity 
of Eqs. (|75jl to give the result in terms of b. This unveils a remarkable similarity of the 
expansions at k = and k = 1. In either case we obtain 

a = a + -^z 2 + a-2-2 4 , & = &o + -^z 2 + &2-2 4 , /? = Po + P2Z 4 , with z = k, k' . (78) 

By abuse of notation we used the same symbols for the expansion coefficients at k = 
and at k — 1. One has to remember that they assume different values depending on the 
phase boundary considered: The constants a , b , /3 obey Eqs. d?H) or (fTBj). respectively. 
They are connected to the value of the effective fermion mass at the phase boundary 
by Mq = a i/so — l//?o or Mq /7 ' ) = a 0y /so/(3o. The coefficients a 2 and 62 become ln(/c')- 
dependent in the case k — ► 1. In either limit they are too complex to write down. However, 
with the integrals 

Jo = Relim f du == — —tanh " 



Ji = Relim / du— — ~ tanh — -, (79) 



we can give the results for /3 2 , 
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Now we can use Eq. or Eq. to calculate ^ryst as a function of k. In order to 
compare the result with ^hom i n Eq. d!ZQ|) [or Eq. (j^B| with k = or k = 1] we have to 
determine the effective mass M for given /3, 7, za Keeping in mind that Eq. connects 
M to the above values of a, s, and (3 at the phase boundary (and only at the phase 
boundary) this is achieved by using c^M^hom = ^ [ or solving F rcn = F™ n = for k = 
or K = 1 to (redundantly) determine the values of a and s in the homogeneous phase 
and thereafter using Eq. (J59*j) ]. We find that the coefficients given in Eq. (J78J) suffice to 
calculate ^ryst — ^hom ^° the fi rs t non-vanishing order, yielding 
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With k 8 = [(j3-ft)/^f and k' 4 = [([3 - (3 ) / (3 2 '] we confirm the behavior of Eq. (J77|). 
[Note that to leading order we have l+41n(4/K') ~ ln(«' ~ 4 ) ~ ln[l/(/3-/3 )] ~ ln(T-T ).] 
We see in Figs. 3, 4 that on a z/ = const, line inside the crystal phase the temperature is 
bounded by its values at the phase boundary: T < T < T . Thus, 

_ agJo r _ 4ao(aoT + fr o7 ) 
^-"512^0/^' * (for 7 >0), (82) 

establishing Cn > 0. Moreover, we deduce > and by Eq. (JHUj) we conclude that 
J < yielding Cj > 0. 

In the massless case 7 = the functional form of Eq. (|77jl is unchanged (see [IH])- In the 
limit 7 -> we have So — > 1 and we observe that Cj converges to the result given in jTH] 
whereas Cn — > AclqTq/tt which is twice the value for 7 = 0. In fact, the chiral limit is 
singular: The k' 4 -coefficients of (3 and of ^ryst — ^Lom diverge in the limit 7-^0 which 
reflects the fact that both expressions are of order k' 2 at 7 = 0. With this singularity we 
cannot expect that Cn behaves smoothly (or even that the functional form of Eq. (|77jl 
prevails). However, the massless and the massive case are still related: In the derivation 
of fi^ 1 ^ the order in k' enters linearly giving rise to the factor of two in the massive case. 

6) The cusp 

The phase boundary sheets I and II meet in a line labelled PtC in Fig. 3. For 7 > 
this line can be determined by taking the limit b — > 0. As discussed in Subsect. 2.4, in 
this limit the scalar potential becomes both homogeneous and ^-independent. For fixed 
v — u c , solving F2 = gives a as a function of b. The effective fermion mass M assumes 
a ^-independent value M c [see Eq. (jfJHj) with s = 1/b 2 + 0(1)], 



where we have introduced the shorthand notation m c and the cusp temperature T c = 1//3 C . 
Because M c is finite, we see that a (b) vanishes in the limit b — > 0, cancelling the 1/b pole 
in S. However, a is not a good expansion parameter at the cusp. Its expansion contains 
negative powers of b [see Eq. ()9(j))] and fails to converge for b — > 0. The situation is 
reversed at 7 = 0. Here b = K in the whole (//, T)-plane and the cusp will be described 
by a = 0. We conclude that at the multicritical point M t = M c | 7= o = connects to 
the massless Fermi gas phase as expected from the phase diagram. We will analyze the 
multicritical point by an expansion at a = in the next subsection. For the rest of this 
subsection we assume 7 > 0. 

By introducing m c and rescaling the integration variable uj —>■ u/b we obtain 



M c = m c T c , 



m c = lim 




(83) 
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By comparison with the homogeneous grand potential of Subsect. 4.3 given by Eq. (J70|) 
one finds (compare Eq. (J72|) for the first equality) 



1 6=0 ' 



d2 <Trrcn 



pren | 



M=M C 1 16=0 ' 

= ~-rrh- (85) 



We see that the minimum equations F± = are equivalent to the conditions \P' = = 
\E r " / = of the translationally invariant calculation [TS|. This proves that the curve b = in 
the new phase diagram coincides with the line of critical points of the old phase diagram. 
We will see next that the opening angle of the cusp is zero (hence the name), just like at 
7 = JH]- Lying inside the crystal phase, the critical first order phase transition sheet of 
the old phase diagram shown in Fig. 2 has then to be tangential to both sheets I and II. 

In order to study the shape of the cusp in detail we expand the minimum equations 
at b = 0. After introducing the physical parameter m c and rescaling u —>■ u/b this is 
straightforward. We found the following behavior 

T-T c = C 1 {v-v c )+C 3/2 {K){v-V c f /2 + C 2 {K){ti- f i c ) 2 + ... , (86) 

n I C T C 1 f°°+ i£ du u _ 2 m c u-v c 
Ci = ; — , i c = -Rehm / cosh . 

4A*c + T/W c — T c 4 e-*0 J-oo+ie \/ UT — 1 2 

The linear coefficient C\ is independent of k which proves that the phase boundaries I 
with k = and II with k = 1 are parallel at the cusp. 

For 7 — » we have b —>■ K in the (//, T)-plane. Hence small values of b are confined to 
smaller and smaller regions about the cusp when 7 goes to zero. We therefore cannot 
expect that the expansion (jBTjl valid for positive 7 converges to the expansion for 7 = 0. 
Comparison with [F3\ actually shows that C\ converges to the chiral value which means 
that the orientation of the cusp changes continuously with 7 — > 0. Moreover, C3/2 vanishes 
in the chiral limit which agrees with the absence of this coefficient for 7 = 0. However, 
C 2 fails to converge to the chiral result. 

The 7-dependence of the cusp can most clearly be seen in the two-dimensional sections 
of Fig. 4. 

7) Multicritical point and Ginzburg- Landau effective action 

The most remarkable of the special points in Table 1, Sect. 4.4, is the multicritical point 
P t at (7 = 0,/i t = 0.608221, T t = 0.318329). At this point the value of a goes to zero and 
we can expand all relevant functions in a. The expansion coefficients are given in terms of 
the polygamma function. This allows us to give an approximate solution of the minimum 
equations F{ = valid in the neighborhood of the critical point for small values of 7. 
Moreover, we will derive a Ginzburg-Landau effective action at the end of this subsection. 

The method to obtain the series in a is identical to the chiral case in [T^]. After a rescaling 
u — > uj/a the expansion in a is obtained from a large 10 series of the algebraic part of the 
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integrand. We define the coefficients c n by a generating function, 



r 

, an 



n=0 



UJ 



We find cq — 1, c\ — (3s — 1 — k 2 )/2. Note that the c n fulfill the recursion relation 
2(n + l)c n+ i = (2n + l)(3s - 1 - /t 2 )c n - 2n(3s 2 - 2s(l + k 2 ) + k 2 )c„_i 



+ (2n - l)(s - l)(s - k 2 )sc 



n-2- 



57) 



Now, we mainly use the residue theorem to obtain the following identities valid in the 
limit A — ► oo 
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(89) 

Expansions of integrands with factor tanh[(ao; — z/)/2] are obtained by differentiation with 
respect to a or v. From Eqs. ()631 ISlj) we derive the following expressions (c_i = 0) 
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These expansions can be used to solve the minimum equations to lowest order in a. The 
result is (one may conveniently use Eq. (|fi5j) for the last identity) 
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The approximation is improved if one corrects TtoT« exp(— 7) / (3 with the above values 
for f3 and 7. If v is close to z/ t then Re ip(2, | + ^) is small and for v > u t the above 
equations give a (6, «, z/)-parameter description of the crystal phase for low 7 near the 
multicritical point. It gives insight into the behavior of the critical line P t C for small 
7. We read off 7 = 0(AT 5//4 ) = 0(A/i 5 / 4 ). A systematic extension to higher orders of 
Re^(2, § + 1^:) is possible. 

We close this subsection with the observation that in the vicinity of the Lifshitz point, 
we can derive a Ginzburg-Landau effective action from Eq. ()9()jl . A similar calculation 
was done at 7 = and found to agree with independent calculations in condensed matter 
physics jU EH] • Here we are interested in a generalization to small but finite values of 
7. Using Eq. (J2U|l one can verify that the coefficients c n+ i — (s — u)c n in Eq. (J9U|) are 
related to the following combinations of spatial averages of powers of S, Eqs. ©, and its 
derivatives, 
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(s e ) + 5(s\sy) + -((s"Y) 



Although the c n and S(x) are different from the corresponding quantities in the chiral 
limit, the fact that Eqs. (|92jl hold enables us to write down the Ginzburg-Landau effective 
action in the simple form 



eff 



eff 



where the effective action at 7 
^ C ff(7 = 0) ■ 



h=o ■ 2-K 
is the same as in Ref. 
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This action can systematically be extended to higher orders in S. The instability with 

respect to crystallization is related to the fact that the "kinetic" (S* 4 + 5" 2 )-term can 
change sign, depending on \i and T. 
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8) Zero temperature limit 



The massive GN model at zero temperature was solved in [T4j . Here we give a short 
account of how to regain the results in our general formalism. 

In the limit T — > both a and v go to infinity. The ratio vja = fi/A, however, tends 
towards a well defined limit which lies in the gap, 

Vs - k 2 < lim - < y/s. (95) 
t^o a 

Moreover, we observe that In 2 cosh[(ac<j — v)j2] — > \{auo — v)/2\ and tanh[(ac<j — v)/2] — > 
sgn(auj — v). From an analytic point of view, this splits the integration contour into 
two sections Keuj > vja and Keu < vja. The orientation of the second section may 
be reversed to compensate for the relative minus sign and then deformed to match the 
first. In addition we obtain a big half circle to account for the singularities at infinity. 
Considering that we are interested in the real part only we may move the left end of the 
contour from vja to \/s, the right end of the gap. This (plus a partial integration for F 2 ) 
transforms the integral representations of \l/ and the Fi to 

,*-Uo = ^-^ + ^(lnA +7 -l)-A 7 (5)-/>-^^- S+M)( ^-^ 



2 2 v ' " ' J^Ts y/W 

limTF ren = -7 + At r dcu^L, 
FH T=0 = 7 + In A - /" du^L, (96) 
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with the limit A — > oo understood. After the substitution uo 2 = z these integrals evaluate 
in a standard way to incomplete elliptic integrals. The result can be extracted from [14j. 



5 Summary and outlook 

The motivation for generalizing the GN model to finite bare fermion masses is twofold. On 
the one hand, from the point of view of strong interaction physics, quarks are massive and 
chiral symmetry is explicitly broken in QCD. On the other hand, a prominent application 
of the GN model in condensed matter physics are conjugated conducting polymers like 
doped polyacetylene. Here, the difference between massless and massive GN models (with 
unbroken or broken discrete chiral symmetry) translates into polymers with degenerate 
or non-degenerate ground states, e.g., trans- and ds-polyacetylene. Since non-degenerate 
polymers prevail in nature, the model with broken symmetry is of considerable interest 
as evidenced by numerous theoretical works on the bipolaron lattice in the 80's and 90's. 

In this paper we have discussed the thermodynamics of the massive GN model in the 
large limit. Our main result is the exact, comprehensive phase diagram of this model 
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as a function of temperature and chemical potential for arbitrary bare fermion mass. 
Previous attempts to solve the same problem have only considered translational invariant 
condensates and produced results which fail theoretical tests. They miss in particular 
the Peierls effect which provides a simple physical explanation for the emergence of a 
crystalline state in one-dimensional fermion systems, namely gap formation at the Fermi 
surface. 

The phase diagram in Fig. 3 brings a series of recent investigations to a close: The T = 
coordinate plane was the subject of Ref. [H] where we addressed the issue of cold, dense 
matter in the massive GN model. The 7 = coordinate plane, i.e., the chiral limit, 
has been investigated first numerically ^7j and later analytically ^H]- The intersection 
of these two planes, the //-axis, was the subject of Ref. [TH]. In these various limits 
integrals which had to be evaluated numerically in the present work simplify. At T = in 
particular, they can be reduced to complete (at 7 = 0) or incomplete (at 7 > 0) elliptic 
integrals. The results in the various limits remain of interest since they give additional 
analytical insights. 

As mentioned before, the GN model is a good example for a fruitful interaction between 
particle and condensed matter physics. At first, semi-classical methods developed for 
baryons in the relativistic case have been taken over to describe solitons and polarons 
in polymer physics. In the following years, condensed matter theory made significant 
progress in the treatment of periodic arrays of polarons. This aspect was overlooked in 
relativistic field theory for a long time, presumably because crystals are not common in 
particle physics. The full phase diagram of the massive model presented here may well 
have applications to real systems in condensed matter physics again. Particle physics and 
condensed matter physics get integrated once one addresses problems of fermionic matter 
and thermodynamics. Incidentally, this trend can also be observed in QCD where interest 
in color superconductivity has led to a similar exchange of ideas [40J. 

Appendix A: Coordinate transformation 

the change of coordinates b, A, k — > (S), (S 2 ),£ in Eq. (JoTlj) has 

: -A(S-U), 

-- -2tA 2 , 

= 0, 

- z + *. 

s 

-- 2A{s - 1 - k 2 + 2u), 
2K 



The Jacobian matrix for 
the following entries, 
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db 

db 
d(S) 
OA 
d(S 2 } 
dA 

m 

OA 
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d(S) A [Z(a - k 2 ) - b(u - k 2 )(s -u)- kH/s] 
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(97) 



The parameters s, t, u are defined in Eq. (ffifj) and Z = Z(b, «) is Jacobi's Zeta function. 
The inverse Jacobian matrix will be denoted by M and its entries are given by 
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-.)(.- rfv • (98) 

In the derivation of Eqs. ()98jl one has to use the relation between t 2 and s given in Eq. ()25jl . 

On the left hand side of Eq. (joTJj) we need the derivatives of ir/3 2 ^/i, Eq. (j52j). with respect 
to the original variables b, A, k. One finds 
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The terms which do not involve an integral arise from differentiating the cutoff A^. In 
the other terms, we differentiate the integrand. When taking the derivative with respect 
to A, we first change integration variables w ^ w/a to avoid having to differentiate the 
ln-factor. Apart from that, we have made repeated use of the structure of the prefactor of 
In (f + e ~ aw + uS ) in the integrand to replace various derivatives by derivatives with respect 
to us. 

We can now invert Eq. (JoTJj) by applying the matrix M of Eq. (jUHj) to the vector (ffilj) and 
adding the contribution from \l/ 2 , Eq. (p)%|) . The result can be cast into the compact form 
given in Eq. (J6U|) as can easily be verified with the help of computer algebra. 



Appendix B: Formulae used in the numerical calcula- 
tions 

In this appendix we transform contour integrals into real integrals of smooth integrands 
with compact support. This step makes the integrals accessible to numerical evaluation. 

The contour of integration of the original integrals is shifted slightly above the real axis, 
with a symmetric cutoff A' at ±00 if needed. The integration variable is u. The integrand 
/ has cuts or singularities on the real line, but no pole at the origin (see Fig. 1 in Subsect. 
3.2). As a function of v it has the symmetry f(—u,i/) = ±f(uj, —u), depending on 
the location of the cuts. We are only interested in the real part of the integral. The 
transformation of the integrals follows the steps (see Fig. 6; some steps may be skipped, 
depending on the integral considered) 



1. Swap gaps with bands: We may swap the cuts with the intervals not cut on the real 
axis without affecting the function in the upper half plane. We thus avoid a branch 
cut at the origin. 

2. Symmetrize the contour: We split the integral into two halves and move the contour 
of one half below the real axis. Due to the branch cuts, we acquire a minus sign 
which we use to invert the orientation of the contour. As a byproduct the imaginary 
parts of the integrals cancel, Re J^ i£ duf = \ (f^ i£ + J^r 1£ )duf. 

3. Cut and paste the contour at the origin: 

rA'+ie , r-A'-ie _ / rO , r- A'-ie\ . / rO . rA'+ie\ 

J-A'+ie ~T JA'-ie — U-A'+ie + JO ) UA'-ie JO )' 

4. Reflect the left path: By symmetry of the contours and a sign change in the measure 
du; we are left with integrating f(u, u)—f(—u, v) along the right path T = A' — is — ► 
— > A' + ie. Using f(—cj, v) = ±/(cj, —u) we have arrived at 

/•A'+ie I r 

Re / f(u>, v) = - / f(u, u) =F (v — > —u), f holomorphic on T. (100) 

J-A'+ie 2 Jr 

5. Isolated poles: Isolated poles are evaluated with the residue theorem. 
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Figure 6: Successive conversion of contour integrals into real integrals with compact 
support. The numbering refers to the corresponding steps explained in the main text. 



6. Shrink the contour: The contour is shrunk to embrace the cuts. 

7. Separate pole contributions: Poles on cuts are subtracted (here, in a symmetric way 
with respect to uj — > — uj) and evaluated separately by elementary integration. Note 
that by the sign change at the cut the subtraction coefficients have opposite signs for 
the arcs above and below the real axis. This is indicated by reversing the orientation 
of the lower arc in Fig. 6. Now that the cuts are free of additional singularities the 
integrand assumes opposite signs above and below the real axis. The value of the 
integral is twice the integral over the upper arc stretching from left to right end of 
the cut, or ending at A' + ie. 

8. Parametrize the cuts: We are left with one or two integrals over (possibly half open) 
cuts. These intervals are parameterized algebraically. For a cut uj e uj\ . . . uj 2 we use 

uj 2 = uj\ + (ul - cj 2 )A 2 with A = — y — and y = . . . 1. (101) 

i + y 

For a half open cut uj G uj 3 . . . A' + ie we take 
uj 2 = uj\ + (uj\ — c<j 2 )A~ 2 with y — 1 . . . 



(A' + ie) 2 -uj 2 -^(A' + ie) 2 -uj 2 



(102) 
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and mi = in case of a single cut. This choice of parameterization lifts the square 
root singularities at the ends of the cut and has the additional advantage to produce 
squares in W, thus simplifying y/W. 

9. Regularize the limit A' — > oo: If necessary we subtract powers of 1/y and evaluate 
these terms in the limit A' —>■ oo by elementary integration. This renders the 
integral finite at y — 0. After swapping the orientation we straighten the contour 
and combine the integrand with a possible second contribution and all constants to 
obtain a single real integral ranging from zero to one. 

10. (optional) Modification for improved numerics: At some instances we have to find 
zeros of functions containing integrals. In this case it is useful to present the function 
as integral over an integrand vanishing at the ends of the integration interval. To 
this end we add the polynomial — f(0)(3y 2 — % + l) — f(l)(3y 2 — 2y) to the integrand 
f{y). Moreover, the integrands may have a lifted singularity which may or may not 
need extra attention depending on the method of numerical integration. In any case 
it is possible to analytically determine the position of the lifted singularity and the 
(finite) value of the integrand at that point. 



In the following we list the integrands f{y), y = ... 1, following the above steps 1 to 9 
(without step 10) for various integrals of 1) the crystal phase, 2) the phase boundary J, 
and 3) the phase boundary II. 

1) Eqs. W!\ \b^\j : The crystal phase 
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tZ+(s-l){s-K 2 ) d u f 



(103) 



2) Eqs. $7$ \7%): Phase boundary I 

sc n A -n tanh ( a _ A _|_ c (A~ 1 y - l)tanhf 

With n — ^= -L - - ^ 

Wltl1 9 "~ 2y 2 s«- 1 (c 2 A- 2 -l) 



a u 

2c 2 
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2s 



In — ^— tanh ( 

2c 2 



n-l 



c = cos 6, s = sin b, we obtain 



'o|«=o : 9o\ a=ac / s , 



A 1 , 

g 2 n H In 2a 



a=a c/s 



, 16s „ 
lim — —Fo : 



Jo 



scd b g 2 \ a=a 



c/s 



S 3 „ C 3 



Ao db J db92 



a=ao c/s, b=bo 



(104) 



5J £gs. |7§ [720: PW boundary II 

sh 2 - n ch"A- 

With k 



y(clrA 



2 A-2 



In 



2 cosh 



/ a 
V2A 



-sh^ln 
2 



2 cosh f 

V2ch 2 



-v) 



ch = cosh 6, sh = sinh 5, 



we obtain 



'o|«-=i 



h 



K=l 



lim In ( — ) Fo : 



C^olo^ach/sh ) 

—d u h 2 In a 

y 



a=a ch/sh 



a sh 



Q=a ch/sh 



Ji : shch d\,d v hi 



Q=ao ch/sh, b=bg 



(105) 
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